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A METHOD FOR EVALUATING SEDIMENTARY BASIN PROPERTIES BY 
NUMERICAL MODELING OF SEDIMENTATION PROCESSES 

CROSS RELATED APPLICATIONS 

[0001] This application claims priority to U.S. Provisional Application No. 
5 ' 60/608,707, which was filed on September 10, 2004. 



FIELD OF THE INVENTION 

[0002] This invention relates generally to the field of geologic modeling and 
reservoir interpretation and characterization. Specifically, the invention is a method 
for constructing a three-dimensional model of the structure and grain-size distribution 
10 of sedimentary rocks that may compose a subsurface hydrocarbon reservoir. 



BACKGROUND OF THE INVENTION 

[0003] A geologic model is a digital representation of the detailed internal 
geometry and rock properties of a subsurface earth volume, such as a petroleum 
reservoir or a sediment-filled basin. In the oil and gas industry, geologic models 

15 provide geologic input to reservoir performance simulations which are used to select 
locations for new wells, estimate hydrocarbon reserves, and plan reservoir- 
development strategies. The spatial distribution of permeability is a key parameter in 
characterizing reservoir performance, and, together with other rock and fluid 
properties, determines the producibility of the reservoir. For sandstone reservoirs, the 

20 spatial distribution of permeability is a function of the grain-size distribution of sands 
which compose the reservoir, the compartmentalization of those sands by barriers of 
finer grained material, and the mineralogy and burial history of the reservoir. 

[0004] Most geologic models built for petroleum applications are in the form of a 
three-dimensional array of model blocks (cells), to which geologic and/or geophysical 
25 properties such as lithology, porosity, acoustic impedance, permeability, and water 
saturation are assigned (such properties will be referred to collectively herein as "rock 
properties"). The entire set of model blocks represents the subsurface earth volume of 
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interest The goal of the geologic-modeling process is to assign rock properties to 
each model block in the geologic model. 

(00051 The geologic modeling process can use many different data types, 
including but not limited to rock-property data obtained from cores, well logs, seismic 

5 data, well test and production data, and structural and stratigraphic surfaces that 
define distinct zones within the model. In general, the resolution or spatial coverage 
of the available data is not adequate to uniquely determine the rock properties in every 
geologic model cell. Other assumptions about the variability in these properties are 
made in order to populate all model cells with reasonable property values. 

10 Geocellular techniques, object-based modeling, and process modeling are three main 
ways to populate the discretized geologic volume with properties. 

[0006] Geocellular models; In the geocellular approach, the relationship 
between properties in nearby cells is specified statistically. Geostatistical estimation 
methods (which may be either deterministic or probabilistic) are used to compute rock 

15 property values within cells. These methods take into account distance, direction, and 
spatial continuity of the rock property being modeled. Deterministic estimation 
methods commonly calculate a minimum-variance estimate of the rock property at 
each block. Probabilistic estimation methods develop distributions of the rock- 
property values and produce a suite of geologic model realizations for the rock 

20 property being modeled, with each realization theoretically being equally probable. 
The spatial continuity of a rock property may be captured by a variogram, a well- 
known technique for quantifying the variability of a rock property as a function of 
separation distance and direction. U.S. Patent Nos. 5,838,634, 6,381,543 and 
6,480,790 cover geocellular modeling methods embodied in processing flows which 

25 include repetitive optimization steps to drive the geologic model toward conformance 
with geologic and geophysical data types such as wells, seismic surveys and 
subsurface fluid production and pressure data. Most commercial modeling software 
packages, including PETREL, GOCAD and STRATAMODEL, contain a wide 
spectrum of geostatistical tools designed to fulfill the requirements of reservoir 

30 geologists and engineers. While these methods can readily accommodate data control 



WO 2006/036389 



PCT/US2005/029884 



-3- 



points such as wells and geophysical constraints such as seismic data, they generally 
do not closely replicate the geologic structures observed in natural systems. 

[0007] Ohiect-based models; In the object-based approach, subsurface reservoir 
volumes are treated as assemblages of geologic objects with pre-defined forms such 

5 as channels and depositional lobes. U.S. Patent No. 6,044,328 discloses one object- 
based modeling scheme that allows geologists and reservoir engineers to select 
geologic objects from an analog library to best match the reservoir being modeled. 
The appropriateness of the analog is judged by the operator of the process based on 
their geologic experience. Most commercial software packages including PETREL, 

10 IRAP-RMS and GOCAD implement objects as volumetric elements that mimic 
channels and lobes using simplified elements based on user deformable shapes such 
as half pipes and ellipses. Other examples of object-based models are the model of 
Mackey and Bridge (1995) and the model of Webb (1994). In their models, the 
depositional objects, such as a river belt in the model of Mackey and Bridge (1995) 

15 and the braided stream network in the model of Webb ( 1 994), are placed sequentially 
on top of each other according to some algorithms. While these models try to mimic 
the real depositional structures, they do not attempt to capture the physics of water 
flow and sediment transport that, over geologic time, determined the rock properties 

t 

at a particular subsurface location. 

20 [0008] Process Models: Process-based models attempt to reproduce subsurface 
stratigraphy by building sedimentary deposits in chronological order relying to 
varying degrees on a model or approximation of the physical processes shaping the 
geology. While process-based models could potentially provide the most accurate 
representation of geologic structures, their application is complicated by the difficulty 

25 in matching their results to data control points or seismic images. Application of 
process models to predicting rock properties generally involves checking the model 
against subsurface data and then rerunning the model with new control variables in an 
iterative process. Although process models are rarely used in current industrial 
practice, U.S. Patent Nos. 5,844,799, 6,205,402 and 6,246,963 describe three such 

30 methods. These methods employ diffusion or rule-based process models to create 
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basin-scale models with limited spatial detail inadequate for reservoir performance 
simulation. 

[0009] A number of academic publications exist on process-based modeling of 
stratigraphy, particularly in shallow water environments. Most of these methods are 

5 not derived from the physics of the depositing fluid flow and sediment transport. 
Rather, transport of sediments is approximated by a diffusion-like equation. Examples 
of models in this class include the model of Bowman and Vail (1999), which is a two- 
dimensional algorithmic model based on a set of empirical rules for a shallow water 
environment, the model of Granjeon and Joseph (1999) which is a three-dimensional 

10 model for simulation of stratigraphy in shallow water environments, the model of 
Kaufman et al. (1991) which has been used in their simulation of sedimentation in 
shallow marine depositional systems, and the model of Ritchie et al. (1999, 2004ab) 
which was used in their work on three-dimensional numerical modeling of deltaic 
depositional sequences. 

15 [00101 Different from the models described above, the SEDSIM model of 
Tetzlaff and Harbaugh (1989) simulates clastic sedimentation in the shallow water 
environment by solving two-dimensional depth-averaged flow equations using the 
Marker-in-cell [Harlow, 1964], or particle-in-cell method [Hockney and Eastwood, 
1981]. Although SEDSIM was originally designed for the shallow water 

20 environment, it has also been applied to simulate a turbidity current by modifying the 
gravitational constant to take into account the relative density of the flow with respect 
to the surrounding medium [Tetzlaff and Harbaugh, 1989]. However, SEDSIM does 
not include the entrainment of water by turbidity currents and the resulting flow 
expansion and impact on sediment deposition. Since SEDSIM uses the Markering- 

25 Cell Method and the flow is represented by many flow elements with constant 
volume, to add the water entrainment relationship would involve significant 
modifications of SEDSIM. 



[00111 The gridding scheme of SEDSIM imposes some limitations. The vertical 
dimension of each three-dimensional cell that contains sediment is fixed. In addition, 
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each cell can only contain a specific single type (or size) of sediment and only the top- 
most cell can participate in the process of erosion. 

(0012J Sediment transport and its effect on the flow are sometimes inaccurately 
modeled in SEDSIM. The reasons include not using experimentally-based erosion 

5 functions, and not allowing, the finer material within the flow to deposit until all the 
coarser material in excess of the transport capacity of the flow has been deposited. In 
addition, SEDSIM does not correct for depositional porosity in modeling of deposit 
elevation so the deposit elevation is underestimated. The effect of spatially and 
temporally varying sediment concentration on the gravitational driving forces is also 

10 not modeled. 

[0013] Syvitski et al. (1999) published a process-based model for simulating the 
movement of sediments onto continental margins and their preservation in the 
stratigraphic record. There are many components in the model of Syvitski et al. 
(1999). They include (1) river mouth process, (2) buoyant river plumes and 
l s hyperpycnal flows, (3) turbidity currents, (4) currents and waves, and (5) debris flows 
generated from slope instabilities. 

[0014] The mathematical model to describe the hyperpycnal flows and turbidity 
currents used by Syvitski et al. (1999) in their model is essentially an improved 
version of the standard Chezy's uniform flow model for a turbidity current [Mulder, 

20 T., J. Syvitski, and K. I. Skene, 1998; Bagnold 1954, 1966; Komar, 1971, 1973, 1977; 
Bowen et al. 1984; Piper and Savoye 1993]. In Chezy's uniform flow model, the sum 
of the gravitational driving force, the frictional resistance to flow, and the internal 
friction of the flow is assumed to be zero. By contrast, in Syvitski et al. (1999), the 
sum of the three forces is (more accurately) assumed to equal the acceleration of the 

25 current. 

[0015] The hyperpycnal flow and turbidity current model used by Syvitski et al. 
(1999) is a large (basin) scale approximation of the flow. As such, it contains many 
assumptions that prevent accurate prediction of the finer-scale sedimentary structures 
that may control the performance of hydrocarbon reservoirs. For example, in their 
30 model, the convection of the flow momentum is not represented. The vast majority of 
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sedimentary bodies that form oil and gas reservoirs are jet and leaf deposits 
[VanWagoner et al. 2003], at which scale the convection of flow momentum plays a 
crucial role in determining the dynamics of the flow and the nature of the bodies. 
Also, the slope of the interface between the turbid and clear water is approximated by 

5 the bed slope in their calculation of the gravitational driving forces. This prevents the 
model from capturing the backwater effect and hydraulic jumps in the flow, which are 
mechanisms in triggering avulsion and filling mini-basins [Beaubouef and Friedman, 
2000]. Without these mechanisms, application of their model in simulating oil and 
gas reservoirs is limited since avulsion is one of the key processes responsible for 

10 reservoir heterogeneities. 

[0016] The hyperpycnal flow and turbidity current model of Svyitski et al. (1999) 
is fundamentally one-dimensional and therefore cannot accurately represent lateral 
variability in rock properties or the physics of the three-dimensional flow that creates 
many characteristic features of deep water sediment stratigraphy as shown by Van 

1 5 Wagoner et al. (2003). Their model describes the confined flow of a turbidity current 
down a one-dimensional conduit. The widths of the conduit have to be specified 
before hand since there is no capability for modeling the initiation and evolution of 
channels and their associated channel deposits, which is another key element in oil 
and gas reservoirs. Finally, the Syvitski et al. (1999) model decouples the interaction 

20 between flow and deposit. 

[0017] In their studies of self-accelerating turbidity currents, Parker et al. (1986) 
derived two sets of mathematic equations to describe the flow of turbidity currents in 
the deep-water environment. The first set of equations is called the "three-equation" 
model. The second set of equations is called the "four-equation" model. Both sets of 

25 equations are depth-averaged equations, and both consist of equations for the 
conservation of flow momentum, flow mass, and sediment. The "four-equation" 
model differs from the "three-equation" model in that it also explicitly takes the 
generation, dissipation, and transport of the turbulent kinetic energy into account. 
Neither model has been combined with the capability to record sedimentary 

30 information about the resulting deposit. 
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[0018] The mathematical formulation of Parker et al. (1986) was originally 
derived for one-dimensional flow with a single sediment grain-size. The "three- 
equation" model has later been extended to 2 dimension and multiple grain sizes by 
Bradford (1996, Ph.D. thesis), Bradford et al. (1997), Imran and Parker (1998) and 
5 Bradford and Katopodes (1999a,b). 

[0019] Both Imran and Parker (1998) and Bradford and Katopodes (1999a,b) 
have applied the extended "three-equation" model of Parker et al. (1986) to study 
turbidity currents and incipient channelization in the deep water environment. In both 
studies, only the bathymetries of deposits in the earliest stages of development [Van 

10 Wagoner et al. 2003] have been reported. Imran and Parker (1998) discussed the 
inability of their model to simulate the longer-term evolution of the bathymetry of the 
deposits, saying "at some point, levee height would approach the thickness of the 
turbidity current itself. The present scheme would fail whenever levees grew to the 
point of nearly completely channelizing the flow, allowing for only small spillover of 

15 the turbidity current." Neither the work of Imran and Parker (1998) nor that of 
Bradford (1997) and Bradford and Katopodes (1999a,b) is related to the simulation of 
the three-dimensional stratigraphy of sedimentary deposits and their associated 
sedimentary bodies and rock properties as shall be described in this invention. 

[0020] Accordingly, there is a need to modify both the "three-equation" model 
20 and the "four-equation" model for the turbidity currents of Parker et al. (1986) to 
simulate the long-term evolution of the sedimentary systems, and to simulate the 
formation of the three-dimensional stratigraphy for deep water deposits. In addition, 
there is a need for a method that honors the shapes and property distributions of 
naturally occurring sedimentary deposits based on geologic data such as seismic and 
25 well data. Such a method may be based on the fundamental laws of physics for water 
and sediment transport and incorporate features from process, object-based and 
geostatistical approaches. Preferably, such a procedure should provide an automated 
optioaso that the optimization process can be performed by a computer, resulting in a 
more accurate model of the subsurface volume of interest, with negligible additional 
30 time and effort. The present invention satisfies this need. 
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SUMMARY OF THE INVENTION 

[00211 A first embodiment of the present invention is a method for simulating the 
formation of sedimentary systems. This embodiment involves, (a) Solving a two- 
dimensional time-dependent map view system of equations for a deep-water turbidity 

5 current, including at least flow momentum, flow height, and suspended sediment 
concentration, and entrainment of overlying clear water, (b) Calculating net sediment 
deposition/erosion at each map view location using the flow properties; (c) Recording 
the net sediment deposition/erosion as it varies with time and spatial location. The 
repetition of this iterative process at each cell of the simulation for each time step 

10 produces a record of sediment deposition history which is a three-dimensional 
geologic model. 

[00221 A second embodiment of the present invention is a method for simulating 
the formation of sedimentary systems. In this embodiment, (a) the initial topography 
of the bottom surface under the simulation region is determined, (b) the locations of 

15 the inlets of the flow that deposited the sediments are determined, (c) a grid is created 
dividing the region of simulation into cells, (d) the flow properties at the inlets are 
estimated, (e) the flow properties in each cell are updated, (f) the net sediment 
deposition/erosion in each cell is calculated, (g) the net sediment deposition/erosion 
for each cell is recorded, (h) steps e-g are repeated until a stopping criterion is 

20 achieved. Optionally, steps d-h can be repeated, adjusting the inlet flow properties 
until the simulated deposit is substantially similar in appearance to an observed 
sedimentary system, thus determining the correct inlet flow properties. 



BRIEF DESCRIPTION OF THE DRAWINGS 

[0023] Figure 1 is a plan view of a fluid flow which is depositing a sedimentary 
25 body, including the flow boundaries; 

[0024] Figure 2 is an elevation view corresponding to plan view figure 1 ; 

[0025] Figure 3 is an elevation view corresponding to plan view figure 1 after 
deposition has occurred; 



t « 
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[00261 Figure 4 is a flow chart of a first embodiment of the invention; 

[0027] Figure 5 is a flow chart of a second embodiment of the invention; 

[00281 Figure 6 is a three-dimensional illustration of an initial bottom surface; 

[00291 Figure 7 is a map view of a simulation in progress indicating the 
5 individual grid cells and the magnitude of the flow velocity; and 

[00301 Figure 8 is a cross-section of a simulated deposit showing bedding and 
grain-size distribution. 



DETAILED DESCRIPTION 

[00311 In the following detailed description, the invention will be described in 
10 connection with its preferred embodiment. However, to the extent that the following 
description is specific to a particular embodiment or a particular use of the invention, 
this is intended to be illustrative only. Accordingly, the invention is not limited to the 
specific embodiment described below, but rather, the invention includes all 
alternatives, modifications, and equivalents falling within the true scope of the 
15 appended claims. 

[00321 The deposition of clastic sedimentary bodies typically begins with a flow 
of sediment-laden water from a confined channel, such as a river mouth or a 
submarine canyon mouth, into an open region, such as a basin. The point where the 
sediment-laden flow enters an open region where deposition occurs is known as the 

20 inlet. Initially such flows expand freely and deposit sediment as the flow decelerates. 
Thereafter, as the deposited sediment grows in height, the deposited sediment begins 
to obstruct the flow field. Eventually, the deposit becomes sufficiently large that the 
flow is diverted around the deposit. This results in a new path and inlet for the flow 
field to an open region beyond or adjacent to the old deposit. The deposition process 

25 then repeats, and a second body in the system is created. In addition, more than one 
such body may be actively built within the system at a time. Overall, the process 
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produces a deposit consisting of stacks of sedimentary bodies, which is the structure 
of hydrocarbon reservoirs. 

[0033] The inventive method enables a physics-based numerical simulation of 
the formation and evolution of three-dimensional sedimentary deposits, their structure 
5 and rock properties. Application of the method can also reduce uncertainties in 
reservoir interpretation and characterization by linking the stratigraphic patterns, 
sedimentary body geometries and rock properties seen in deposits with the associated 
depositional processes and the corresponding geologic controls in the deep-water 
environment. 

10 [00341 The inventive method is based on the physical laws of fluid flow and 
sediment transport. It is capable of simulating the large scale and long term evolution 
of depositional systems, including the development of the associated sedimentary 
structures and rock properties, with high spatial resolution in three-dimensions. The 
method, depending on the embodiment, may involve (1) using a physics-based 

15 numerical simulation for the calculation of fluid flow and sediment transport in the 
deep water environment with given initial and boundary conditions; (2) coupling the 
fluid flow and sediment transport model with an appropriate erosion and deposition 
model to simulate the erosion and deposition of sediments of multiple sizes and 
properties everywhere in the system; (3) using an appropriate three-dimensional 

20 griding scheme to record the erosion and the deposition of sediment everywhere in the 
system, as well as the changing bathymetry of the system; (4) asserting external 
controls on the system by changing the boundary conditions dynamically. The 
combination of these features leads to a method for constructing a geologic model that 
specifies the grain-size distribution throughout a subsurface volume, guided by a 

25 seismic image of the subsurface volume. This method, in its various embodiments, is 
the subject of the present invention. 

[00351 Figures 1, 2, and 3 depict the assumptions and parameters used in the 
present invention. Figures 4 and 5 are flow chart diagrams of two embodiments of the 
present invention. Table 1 provides a detailed list of the variables used in the present 
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method with references to the figures, when applicable. The figures are described in 
greater detail below. 

■ 

[0036] Figure 1 depicts a plan view of fluid flow 10 with flow boundaries 12 
and 14. Inlet 15 for fluid flow 10 is centered, for convenience, at the origin of the x 
5 and y-axes, and flow emitted from the inlet moves initially in the positive x direction. 
At the inlet 15, the flow boundary has an initial width 17 and a half-width 8, and 
expands in the positive x direction. Also depicted is an outline of the deposit 16 

* 

formed by the flow. 

[00371 Figure 2 depicts an elevation view corresponding to plan view Figure 1 
10 along the x-axis. The fluid inside the flow boundaries 12 and 14 of Figure 1 is 

comprised of two layers. Figure 2 illustrates the two layers of fluid as a clear layer 20 

above a sediment-laden layer 28 and bounded by the sea surface 21 and sea floor 22. 

The sediment-laden layer 28 is also referred to herein as the turbid water layer. The 

elevation of the bottom before the deposition process occurs is 24. The height of the 
15 sediment-laden water layer at inlet 15 is 26. The height of the sediment-laden water 

layer 28 may vary based on location, as illustrated by height 27 located at a different 

point further along the x-axis in Figure 2. 

[0038] Figure 3 is an elevation view corresponding to the plan view of Figure 1 
along the x-axis after deposition has occurred. As in Figure 2, the fluid flow in Figure 
20 3 is depicted as being comprised of clear layer 20 above a sediment-laden layer 28. 
The elevation of the bottom after deposition is 30. This elevation consists of the 
elevation of the original bottom 24 together with the thickness of the newly deposited 
sediment layer 36. 

[0039] Figure 4 depicts a flow-chart diagram of one embodiment of the present 
25 invention and Figure 5 depicts a flow-chart diagram of a second embodiment of the 
present invention. 
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MODEL EQUATIONS 

[0040] In the model, the mathematical equations that govern the flow of the fluid 
and the transport of the sediments are depth-averaged time-dependent flow equations 
based on the previous work of Parker et al. (1986), Bradford (1996, Ph.D. thesis), 

5 Bradford et al. (1997), Imran and Parker (1998) and Bradford and Katopodes 
(1999a,b). However, the inventors have made some additions and changes to the 
model which enable the numerical simulation to better capture various features of the 
depositional process, such as avulsion and channelization, as well as to better match 
flume tank experiment results. The mathematical model used in the inventive method 

10 and specific details for numerical solution of its equations are described below. 

1 ) Momentum. Sediment Fluid, and Turbulent Ener gy Balance equations 

An equation for the balance of the x component momentum is: 

dt dx dy 2 dx dx dx [l] 

- c d i]u 2 x +u 2 y u x + r x - u x S w sgn(£ w - s w ) 

The equation for the balance of the y component momentum is: 



15 



dt dx dy 2 ' 6 dy 1 dy dy [2] 

An equation for the conservation of the turbid water layer: 

dh duji 3* , +r [3] 
dt dx dy 



An equation for the conservation of the sediments with multiple grain-sizes in the 
flow is: 



20 
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dhC, + dhuA + m £ ^ for , m 1>2 „ [4] 

dt dx dy 

[0041] Optionally, an additional equation describing the balance of turbulent 
kinetic energy can be included. Here the original one-dimensional single grain-size 
formulation of Parker (1986) has been extended to a two-dimensional multiple grain- 
5 size form. The modified formulation for K is: 



10 



Fit Br dv 2 M rjj 



a/ a* ay 

--a,RgChJ^Tule w -^a t Rgh^ l v t ,(E,-r 0 C l ) + t K -Kd w 8sa(.S w '8 w ) 

[0042] In the above equations, « x and u y are the x - and y - components of the 
depth-averaged flow velocity respectively, h is the depth of the turbid flow, t) is the 
bottom elevation, C, is the volumetric concentration of the sediments in the i th 

grain-size bin andC = £C, is the total sediment concentration. In equations [1] and 



[2], g is the gravitational constant, and R is the submerged specific weight for the 
sediments where R = P ' ~ Pw and p, and p w are the sediment and water density, 

Py, 

respectively. In equation [5], «.is the shear velocity, where ul =a K K , and a K is a 
model parameter related to the turbulence. The remaining terms in equations [1] - [4] 

15 and [5] may be further described below, but are identified briefly here. In equations 
[l]-[2], a, is the stratification parameter which characterizes the vertical variation of 
sediment concentration within the flow. In equation [5], P K is another parameter of 
the turbulence. In equations [l]-[2], c d is the friction coefficient producing a 
velocity-dependent drag on the flow. In equation [3], e w is the entrainment function 

20 and ^is the detrainment function. The entrainment function characterizes the rate at 
which stationary clear water above is entrained into the moving turbid water below, 
becoming part of the flow. The detrainment function characterizes the reduction of 
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the flow height due to the setting of the sediment from the top most part of the 
turbidity flow (Toniolio, 2002). In equation [4], E, and D pi are the erosion and 
deposition functions that characterize the rate of erosion of sediment in the i th grain- 
size bin from the bottom into the flow and the rate of deposition of sediment in the 
5 j th grain-size bin to the bottom from the flow, respectively. In equations [1], [2], and 
[4], Tx, Ty, t w r c and r K are optional terms describing the transport due to turbulent 
diffusion. The optional term sgn( ), which is shown in equations [1], [2], [3], and [5] 
is defined as sgn(a) = 1 for a > 0 and sgn(a) = 0 for a < 0. 

(2) The stratification coefficient ou 

10 

[0043] In one embodiment, the inventive method modifies the previous equations 
of Parker et al. (1986), Bradford (1996, Ph.D. thesis), Bradford et al. (1997), Imran 
and Parker (1998) and Bradford and Katopodes (1999a,b) by introducing a 
stratification parameter a t in the momentum conservation equations of [1] and [2]. 
15 Sediment tends to be dominantly carried in the lower portion of turbidity flows, not 
distributed uniformly through the thickness of the flow as the earlier work had 
assumed. Assuming a stratification parameter value less than 1 accounts for the fact 
that the center of mass of the flow is generally below half the height of the flow. 

[0044] When the stratification effect has been taken into account by using a value 
20 of a, < 1 , the component of the gravitational force associated with the gradients of 
both flow height and sediment concentration are reduced. As a result, the component 
of the gravitational force associated with the gradient of the bottom elevation becomes 
more significant. Therefore, the stratified flow will have stronger interactions with 
the underlying topography and bathymetry. The inventors discovered that without 
25 using this correction, it becomes difficult to accurately model flow avulsion and the 
model becomes more susceptible to numerical instabilities. When the effect of 
stratification has been taken into account, the model results match more closely the 
results from flume experiments. 



* • 
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[00451 The values of a s to be used in simulations can be determined in a number 
of ways. An approximate value of a, can be selected, typically 2/3. Alternately, the 
value of a, can be adjusted by the user such that the overall flow and sediment 
transport properties and their associated overall sediment deposition (or erosion) 
patterns best match the available experimental or field data sets. 

[00461 While a constant value of a, may be used in the model, a, can also be 
made to depend on other flow parameters such as the local Froude number Fr , as, for 
example, to use one value ofa, (,) for a, when the Froude number of the flow F r is 
< 1 and the other value a, (2) for a s when F r > 1 . Both of these values can be user 
10 input parameters and can be adjusted to make the simulation correspond to observed 
data from the natural system. 

(00471 Another method to determine values for a, is based on a modification of 
the "top-hat" or "slab" assumption for the vertical sediment concentration profile 
[Turner (1973) and Pantin (1979)]. m Parker's original derivation [Parker, 1986], the 
1 5 following "top-hat" or "slab" assumption was used. 

flfor0<^<l rfil 
fc<H>«ff.(H)-fcr<*>-| O for„>l ' 1 J 

In above equation, g a (rj) , $M and g k (v) are the vertical profiles for the flow 
velocity, sediment concentration and turbulent kinetic energy. In the same equation, 

tj = - is the rescaled (by flow depth h ) dimensionless z-coordinate. Different from 
h 

20 Parker's original model, we use: 



fl forO<77^1 
f.fo)-ffr(l>-| O for^>l 



[7] 

/ c ( 7 )for0<^<l 
0 for77 >1 
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where fArj) is some normalized bottom-weighted function that captures nature of the 



sediment stratification. The normalization is, 



[8] 



The vertical profiles we proposed in equations [7] and [8] are consistent with the 
5 constraints set by the original derivation of the governing equations as shown in 
equations [1] to [4]. These constrains are, 



1 



10 In this method, a s can then be calculated as, 



[9] 



For non-stratified flow, a, = 1 . 



15 [0048J The choice of f c (g) includes, but is not limited to, the following two 
examples: 



(a)/ c («)-(r+D0-rt' I 10 ' 



20 
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where in above equations, both yand £ 0 are model parameters that can be used to 
adjust the degree of stratification. It will be shown below however, that neither the 
explicit form of f c (g) , nor the specific choice of y and ^are utilized in the final 
form of our simulation models. They will be encapsulated into a single parameter a, , 
which we will term the geometric factor. By choosing appropriate values of a, , 
varying degrees of stratification effects can be modeled. 



(3) Turbulence parameter s o^, Bkt 

10 [0049] A difference between one preferred embodiment of the inventive method 
and Parker's original formulation lies in the treatment of the parameter/?*, in equation 
[5], In Parker (1985), /?*has been related to a characteristic friction coefficient c d 
based on the assumption that at steady state, the friction coefficient c d is the same 
everywhere and has a value equal to c d . In our model however, we do not make the 

15 same assumption and consequently,/?, is specified directly and is treated as a model 
parameter. The range of fi K is generally between 0.0001 and 1 or even 0.005 and 0.5, 
with a typical value being 0.05. 

[0050] The turbulence parameter a K is supplied by the user. The range of a K is 
generally between 0.001 and 1 or 0.0001 and 1, with a typical value being 0.01. 

20 f <tt The friction coefficient Ca. 

[0051] The friction coefficient c d can be treated as a model input parameter. The 
most common range of c d in both experimental and field observations is from 0.001 
to 0. 1 or 0.002 to 0.05, with field experiments using lower values, typically 0.004, and 
lab experiments using higher values, typically 0.02. The friction coefficient can also 
25 be allowed to vary spatially, where the function c d (x,y,z) may be defined by the user 
or may be linked to local flow properties and the associated turbulence. For example, 
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c d may be defined as a function of local Reynolds number. Alternatively, Parker 
(1986) has related c d to the local turbulent kinetic energy K and flow velocity via, 



[12] 



2 

y 



where a K is the turbulent parameter discussed previously. In order to use this last 
5 relationship, the creation, transport, and dissipation of turbulent kinetic energy must 
be modeled explicitly, as in equation [5]. 

(5) Entrapment and detrainment 

[00521 The entrainment function describes the entraining of clear water by the 
flowing turbid water. There are many different forms of the entrainment function. 
1 0 Any of them can be used in this model. The entrainment function used by Parker et al. 
(1986) is 



£... = 



0.00153 



w 0.0204 + R, ' 



[13] 



where fl,is the Richardson number and equals to the inverse of the square root of the 
well known Froude number F r , namely 



15 



F 



[14] 



and 



1^ 1^ y 



[15] 



The same entrainment function has also been used in the current model 
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[0053] Opposite to the entrainment of the clear water into the turbidity flow, 
settling of the sediment from the topmost part of the flow detrains clear water from 
the turbidity current and return it back to the surrounding environment. The concept 
of detrainment was first proposed by Toniolo et al. (2002). The original formulation is 
5 for a single grain-size and relates the rate of detrainment to the grain settling velocity 
S = v . The formulation we utilized is based on the original ideas of Toniolo et al. 
with extensions to make them applicable to turbidity currents carrying sediments of 
multiple grain-sizes, where we propose 

*„ = v,(2>'), t 16 l 

10 and D* is the effective grain-size that characterizes the overall settling interface of the 
turbidity flow. The actual value of D* could range from the minimum grain-size to 
the geometric mean grain-size of the sediment present in the flow. An example of the 

choices for D* is 



D' =D 



10 



[17] 



15 where D l0 is the diameter of the 10/A percentile in the distribution. 
(6) Deposition rate and settling velocity 

The deposition rate is typically the deposition rate of sediment in still water. 

A = r 0 C t v Ml 



[18] 



where r 0 is a model coefficient that relates the bulk sediment concentration to the near 
20 bed concentration, and v„ is the settling velocity of the sediments in grain-size bin i . 

[0054] The settling velocity v,(D) for a sediment grain with diameter Dean be 
specified in a number of different ways. Four methods are described below as 
examples. 
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[0055] The first method is to use the Dietrich's settling velocity curve [Dietrich 
(1982)] for v,(D). In addition to the grain-diameter D, Dietrich has also provided 
corrections of v,(D) for the geometric and material properties of the grain, such as 
the aspect ratio of the grain, the surface roughness of the grain and the flatness of the 
5 grain. Although Dietrich's relationship assumes no interactions among the particles 
in suspension, it is a good approximation for sediment in most natural turbidity 
currents where the concentrations are typically low. 

[0056] The second method is to base the value of v,(D) mostly on Dietrich's 
settling velocity curve v, lDkMi {D) , but apply corrections to the settling velocities of 
10 the finest material in order to take flocculation into account. One of the simplest 
method to take the effect of flocculation into account is to modify the settling velocity 
as follows: 



v r rieh \D),iiD>D M [19] 

where Djioc is the effective diameter of the settling floe. 

15 [0057] The third method is to use Dietrich's settling velocity curve as the base, 
and apply corrections to take the effect of the hindered settling into account. 
Hindered settling occurs at high sediment concentrations where interactions among 
suspended particles produce settling velocities significantly different from those of 
free-settling particles. While most of the sedimentary rocks are believed to be 

20 deposited by flows with relatively low concentrations where the fall of the sediments 
from the flow to the underlying bed can be relatively closely approximated by free 
settling particles, an increasing body of evidence has shown that there are many 
sedimentary rocks deposited by high concentration flows in which the fall of the 
particles from the flow can only be described accurately by including the effect of the 

25 hindered settling. In one embodiment, the effect of the hindered settling can be taken 
into account by applying a hindering function, which corrects the particle settling 
velocity from that of the free settling fall velocity, 
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[20] 



» 

[0058] One example of such a hindering function is the polydisperse hindering 
function of Patwardhan and Tien (1985). The hindering function of Patwardhan and 
Tien (1985) is based on a one-dimensional fully turbulent sedimentation model. In 
that function, the hindering effects are a function of the sorting of the grain-size 
mixture as well as bulk concentration. In poorly sorted suspensions (SD > 1 phi unit) 
hindering may occur at concentrations lower than 0.5%. In well-sorted suspensions 
(SD < 1 phi unit) hindering is negligible at 5% and higher depending on the sorting. 

[0059] The hindering function w w „ rf of Patwardhan and Tien (1985) is, 



10 



1- 



1 



A,[(i-<*>r' /3 -i 

A 



i+ 



6,-2 



[21] 



In above equation, O is the total porosity in the flow, D m is the geometric mean grain- 
size suspended in the flow, 2} is the grain-size of the sedimentary particles in the i th 
bin, and 0, is a function of the particle Reynolds number 



0, = 4.35*,, 
0, = 4.45J?„ 



-0.03 



0.2 < R pl < 1 
1 < R pl < 500 



[22] 



15 The particle Reynolds number R pi for the grains in the / th size bin is defined as 



[23] 



where v is the kinematic viscosity of the water 
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(7\ Erosion rate. 

[00601 Examples of the erosion functions that can be used in this model are given 
here. The most commonly used erosion function is the Garcia-Parker erosion 
function [Garcia and Parker, 1991]. In this function, the rate of entrainment of 
5 sediment of grain-size bin i into the flow from the bed is 



10 



= a ' Z ' v,,G, [24] 



a. „ 5 



1 + ^-Z: 



where Z, is a function defined as 



x0.2 

D 1 [25] 



in which 



Ir_, 06 if i?„ f > 2.36 [26] 



f( * Rp,) {o.586i?„'- 23 if R pi < 2.36 



and 



X = \- 0.288a . t 27 ] 



[0061] In equation [24], G t is the volumetric percentage of the sediments of 
grain-size bin i in the surface layer, a, is a constant and typically has a value of 

£ 

15 1 3 x 10" 7 , e is the maximum value of the dimensionless erosion rate — and it sets 

v ti G, 

the upper limit of the erosion function. In equation [25], £> J0 is the diameter of the 
sediment grain in the 50th percentile in the distribution. In equation [27], a is the 
standard deviation of the grain-size distribution in the logarithmic "phi" units familiar 
to geologists. 
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[00621 An immediate variation of the erosion function is to calculate Z, as 



u. 



[28] 



si 



rather than using the function shown in equation [25] and to set X = 1 instead of using 
equation [27]. This simplified erosion function is found to fit some of flume 
experiments better than the first erosion function. 

[00631 A third example of an erosion function that could be used is from 
Akiyama and Fukushima (1985). In this erosion function, the rate of entrainment of 
sediment of grain-size bin i into the flow is: 



0.3v„G,, Z,>Z m 

o, z'<z c 



Z c <Z i <,Z m 



[29] 



10 



where Z c = 5 and Z m = 13.2 . This erosion function was used in the study of Parker 



et al. (1986). 



(S) Turbulent diffusion. t» t v , t c 



[0064] The following terms are added to the right side of equations [1 ], [2] and [4] 
respectively to model the spreading effect caused by turbulent diffusion. Turbulent 
15 diffusion is the time averaged transport of mass and momentum by turbulent eddies. 
These three terms are defined as follows: 

(a). The term for the turbulent diffusion of the x-component of momentum 



Tv = 



dx 



dx J dy 



dy J 



[30a] 



(b). The term for the turbulent diffusion of the y-component of momentum 
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d 


r d(u y hy 


d 

+ — 


d(u y h) 




v — - — 


dy 


v, — - — 



[30b] 



(c). The term for the turbulent diffusion of the sediment concentration for each 
grain-size bin i 



Tci = 



dx 



d(hc f ) 
dx J dy 



+ 



. d(hc,) 
' dy 



[31a] 



(d). The term for the turbulent diffusion of the water 



10 



dx 



dh 

/' dxV dy 



dh 



[31b] 



(e). The term for the turbulent diffusion of the turbulent kinetic energy 



3Kk] + d 



dx J dy\_ dy _ 




dKh 



[32] 



15 



In the above terms, the turbulent diffusion coefficient v, can be related to the von 
Karman constant k, shear velocity «.of the flow, an empirical correction factor, A, 
which ranges from 10 to 100, and is typically 30, and the flow depth h by the 
following equation, 



20 



v, = A—u.h. 
6 



[33] 



While the above three terms are not new for shallow-water equations, an embodiment 
of the inventive method is the first time that they have been applied to turbidity 
currents in association of the set of the governing equations as shown in equations [1- 

4]. 
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NTTMPRTCAT, SOLUTI ON OF THE MQ DF.T. EQUATIONS 
(tt Gridding 

[0065] To solve the equations numerically, the region of simulation is divided 
into cells. One gridding technique which is well suited for these simulations is 

5 described in co-pending U.S. Provisional Patent Application No. 60/584,617 titled 
"Method for geological modeling through hydrodynamics-based gridding (Hydro- 
Grids)." In Hydro-Grids, horizontal grid lines in the deposit model correspond to 
surfaces of constant time during the simulation. Vertical grid lines are also created 
based on rates of structural variability in available geologic data. The sedimentary 

10 properties of the deposit are then efficiently represented by these hydrodynamics- 
based grid cells. 

(2\ Solver 

[0066J The numerical solution of the model equations is complicated because of 
numerical instabilities. The choice of solver appears to be used in avoiding these 

15 instabilities. The solver is a method to obtain numerical solutions for the partial 
differential equations in a discrete domain. Imran, et al. (1998) reported that their 
model crashes consistently at the beginning of flow avulsions. Bradford and 
Katopodes (1999) used Roe's solver [Roe, 1981] in their studies. However, none of 
their simulations have passed the stage of incipient avulsion. The inventors' own test 

20 of Roe's solver has shown similar numerical instabilities to those reported by Imran et 

■ 

al. (1998) developing as the simulation approaches the stage when avulsion would 
occur. 

[0067] An innovative aspect of the preferred embodiment of the inventive 
method is using an HLL solver to solve equations [1] to [5]. Different from Roe's 
25 solver, the HLL solver obtains an approximate solution for the inter-cell numerical 
flux directly by assuming a particular wave configuration for the solution and 
calculating the wave speed [Harten et al., 1983]. Although the HLL solver has been 
used to solve shallow water equations before, it has never been used for the depth- 
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averaged equations for turbidity currents as shown in equations [1] to [5]. The HLL 
solver has proven stable even when running long enough to generate highly 
complicated stratigraphy. 

[00681 The solver steps the solution of the equations forward in time. At each 
5 time step, the flow properties in each cell are updated, the deposition rates are 
calculated, and the bottom is modified to reflect the deposition. The solver completes 
when some stopping condition is reached. Typically the stopping condition is the 
geologic time duration to be simulated. Alternatively, the stopping condition could be 
based on the deposited sediments reaching some thickness. 



10 (3) Boundary Conditions 

[0069] Boundary conditions which can be implemented in one embodiment of 
the inventive method include: standard open boundary conditions, wall boundary 
conditions, and dry nodes. In addition to cells which are defined to be dry as a 
boundary condition, some cells that were previously covered by flows (either turbidity 
15 currents in the deep-water model, or just water flow in the shallow water model) 
could emerge from the flow and become "dry" cells due to a variety of reasons, for 
example, upstream avulsion of the flow. Emerging "dry" nodes often represent 
sources of numerical instabilities and create inaccuracies in the calculation. The 
situation becomes especially severe in regions of high bed topography gradient. 

20 [0070] The algorithm we propose to solve the problem associated with "dry" 
nodes is described below. Without loss of generality , assume two neighboring cells A 
and B whose bed elevation are tj a and rj B respectively. The flow height on top of 
these two cells are h A and h B respectively. 



If h A > 0 and Tj B +h B <r} A , the flux between the two cells is calculated using 
25 following algorithm: 
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(1 ) Calculate the flux F AB in Cell A at the edge that bonds Cell B by 
temporarily setting the flow depth h B in Cell B to zero. 



(2) Calculate the flux F B in Cell B at the edge that bonds Cell A by 
temporarily setting tj a to oo . 



(3) The out flux in Cell A at the edge that bonds Cell B is F AB . 



(4) The out flux in Cell B at the edge that bonds Cell A is then F B - F AB . 



When h A = 0 and tj b + h B < r\ A , the flux between the two cells is calculated using 
following algorithm: 



(1) Calculate the flux F B in Cell B at the edge that bonds Cell A by 
10 temporarily setting ij A to oo. 



(2) The out flux in Cells A and B is zero and F B , respectively. 



(41 Speed enhancement bv reducing number of grain-size bins 



[00711 In computer simulation, the continuous distribution of sediment grain-size 
is represented by a finite number of grain-size bins. Each bin has an associated 
1 5 nominal grain-diameter, typically the geometric mean of the upper and lower limits of 
the grain-size range contained in the bin. All sediment within the bin is assumed to 
erode and deposit like sediment having this nominal grain-diameter. 

[00721 For computational speed, it is desirable to limit the number of grain-size 
bins. However, if the number of bins is too small, the sediment within a bin may be 
20 dominantly larger or smaller than the nominal diameter. In this case, the net 
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deposition rate for the bin will be miscalculated, and the size and shape of the 
simulated deposit will be compromised. 

[0073] This problem can be mitigated by dynamically adjusting the nominal 
grain-diameter for each bin so that the net deposition rate is correctly predicted. With 

5 this novel correction, the number of bins utilized for accurate simulation can be 
reduced. In one embodiment, the appropriate nominal grain-diameters for each bin 
are found assuming that a) the complete grain-size distribution in the flow is log- 
normal in grain-diameter and b) the net deposition rate within each bin is exponential 
in grain-diameter. Other methods for computing appropriate nominal grain-diameters 

10 based on interpolating the grain-size distribution and the net deposition rate within 
bins are also within the scope of this invention. 

[0074] In this discussion, grain-diameter x will be given on the logarithmic phi 
scale familiar to sedimentologists. The log-normal grain-size distribution in the flow 
can then be written as, 



15 p(x) = -7=!=-exp| 

V2^(7 



(*-/0 2 
2a 2 



[35] 



where p(x) is the differential volume fraction of grains having diameter x, n is the 
mean grain-diameter on the phi scale, and a is the standard deviation of grain- 
diameter on the phi scale. The mean and standard deviation of this distribution at a 
given location and time step are determined by curve fitting to the fractions of the 
20 total grain volume within each size bin. 

[0075] Within the ith grain-size bin, the net deposition y can be written as a 
linear function of the phi scale grain-diameter, 



y = a,x + b,, 



[36] 



where a, and b ( are computed to match the net deposition rate for the largest and 
25 smallest grains included in the ith grain-size bin. 
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[0076] Using equation [35] and [36], the total net deposition Y t for the ith grain- 
size bin can be expressed as follows, 



^\a i x + b l )p(x)dx 



£p(x)dx 




exp 


f 1 / 

L 2a L J 


-exp 


1 / \2 
{ 2cr z J 


erf 











[37] 



a, +6,. 



[0077] From equation [37], the effective grain-size X (on the phi scale) which 
will produce the same net deposition as the continuous grain-size distribution, can be 
estimated using the following formula: 



• » 

WO 2006/036389 



30 



PCT/US2005/029884 



EMBODIMENTS 

[0078] A first embodiment of the invention will now be described. With 
reference to Figure 4, this embodiment is a method for simulating the formation of 
sedimentary systems. As illustrated in Figure 4, this embodiment involves solving a 

5 two-dimensional time-dependent map view system of equations for the flow of 
sediment-laden water (step 401), calculating net sediment deposition at each map 
view location using the flow properties of the two-dimensional map view system of 
equations (step 402), recording the net sediment deposition as it varies with time and 
spatial location (step 403), and repeating, if necessary, steps 401-403 for each time 

10 step until the time duration of the simulation is complete. The individual steps will be 
described in greater detail in the following paragraphs. 

[00791 In one embodiment, step 401 involves determining the properties of the 
turbidity flow by solving a two-dimensional time-dependent map view system of 
equations, where the flow properties described by the equations comprise flow 

15 momentum, flow height, suspended sediment concentration, and entrainment of 
overlying clear water. This solution process expresses equations [l]-[4] in discrete 
form where the value of the flow parameters in each cell are updated according to the 
method specified for the HLL solver [Harten et al. 1983]. To use the HLL solver, 
equations [l]-[4] are first written in a conservative form, where the primary variables 

20 such flow height, flow velocity and the suspended sediment concentration are 
transformed to a set of conservative variables defined as follows: 

[0080] U x = u x h . 

[0081] U 2 = u y h . 

[0082] U 3 =h. 

25 [0083] C/ 4 , = C,/i where i = 1,2,...,«. 

[0084] Once the values of the these conservative variables are determined, the 
primary flow variables can be easily calculated using the following equations: 
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[00851 u x =j± 



100861 u y =jf- 



[00871 h = U 3 . 



[00881 C t = — - where i = 1 ,2,3,..., n . 

5 [00891 At each time step, the values of the conservative variables U ]t U 2 ,U 3 
and C/. are evaluated for each cell. In the HLL solver, the conservative variables are 
assumed to be piecewise constant in each cell. The conservative variables are also 
assumed to be discontinuous across all cell edges that separate two cells, and the 
discontinuities are treated as shock fronts. The wave speed associated with the 

10 propagation of the shock fronts across all the cell faces are then calculated. Once the 
wave propagation speeds are obtained, fluxes across that cell face are then calculated 
for each conservative variable. The values of the conservative variables in each cell 
are then updated, and the new values for the flow parameters in each cell are obtained 
from these newly updated conservative variables. 

15 [00901 Step 402 involves calculating net sediment deposition at the current time 
step for at least one map view location using the flow properties updated in step 401. 
The properties of the net sediment deposition that can be calculated include but are 
not limited to the net deposition rate for each grain size bin, the total thickness of 
sediment deposited during the time step, the new bottom elevation at the cell, bedding 

20 type within the cell (as inferred using flow velocity), and any combination thereof. In 
one embodiment this step may be performed using equations [18] and [24]. 

[00911 Step 403 involves recording the net sediment deposition at the current 
time step for at least one map view location. The net sediment deposition, as it varies 
with time and spatial location, may be updated and recorded in a grid of model cells 
25 during the simulation. One gridding method suitable for recording the net sediment 
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deposition at it varies with time and spatial location is described in co-pending U.S. 
Provisional Patent Application No. 60/584,617 (Hydro-Grids). 

[0092] In an optional additional step, Steps 401 through 403 are repeated for each 
desired time step until the time duration of the simulation is complete. Furthermore, 
5 the simulated sedimentary deposit may be checked against available geologic data 
such as seismic, outcrop, core, or well log data. If the simulated deposit does not 
match the geologic data, boundary conditions (including inlet flow properties and 
initial bottom topography) may be adjusted and the steps repeated in an iterative 
manner. 

10 [0093] A second embodiment will now be described. With reference to Figure 5, 
this embodiment is a method for simulating the formation of sedimentary systems. As 
illustrated in Figure 5, the topography of the surface underlying the simulation region 
is determined (step 501), the location of the inlet of the flow that deposited the 
sediments is determined (step 502), a grid is created dividing the region of simulation 

15 into cells (step 503), the flow properties at the inlet are estimated (step 504), a flow 
simulation is performed by updating the flow properties in each grid cell based on the 
estimated flow properties at the inlet (step 505), the net sediment deposition in each 
cell is calculated based on the flow properties (step 506), the net sediment deposition 
for each cell is recorded (step 507), steps 505-507 are repeated until a stopping 

20 criterion is achieved (step 508). The individual steps will be described in greater detail 
in the following paragraphs. 

[0094] Step 501 involves determining the topography of the surface underlying 
the simulation region. In one embodiment this is typically accomplished by 
identifying a stratigraphic surface in a three-dimensional seismic data volume. This 
25 can be accomplished by methods familiar to persons of ordinary skill in the art, who 
will also recognize other methods for identifying or inferring the geometry of a 
stratigraphic surface. Such other methods include, but are not limited to interpretation 
of two-dimensional seismic lines, other remote imaging techniques, correlating well 
logs, and spatially correlating outcrop observations. 
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[0095] Step 502 involves determining the location of the inlet of the flow that 
deposited the sediments above the bottom surface. This is typically accomplished 
based on the interpreted paleoflow direction determined by a seismic interpreter 
according to methods familiar to persons of ordinary skill in the art. 

10096] Step 503. A grid is created by dividing the region of simulation into cells. 
The preferred embodiment of the gridding method is described in co-pending U.S. 
Provisional Patent Application No. 60/584,617. 

(00971 Step 504. The inlet properties are estimated. The inlet properties 
characterized in this step comprise the flow velocity, the flow height, the volume 
10 fraction of the flow composed of grains within each grain-size range, and the inlet 
width. Optionally, corrections to the bottom surface can be estimated at this step as 
well. Such corrections could include overall slope adjustments to the surface to 
restore the paleoslope at the time of deposition. 

[00981 In co-pending U.S. Provisional Patent Application No. 60/466,655, 
15 assumptions are made about relationships between the inlet flow properties. In one 
embodiment, flow velocity becomes the only independent variable in defining the 
local inlet properties. Thus, in that embodiment, characterizing the inlet properties 
can be accomplished by simply estimating the inlet flow velocity. 

[0099] Step 505 involves simulating deposition by updating the flow properties 
20 in each grid cell based on the estimated flow properties at the inlet. Equations [l]-[4] 
are written in discrete form. Therefore, at each time step the flow parameter values in 
each cell are updated according to the method specified for the HLL solver as 
described previously [Harten et al., 1983] 

[0100] Step 506 involves calculating the net sediment deposition for each cell 
25 from the flow parameters. Typically equations [18] and [24] are used to calculate the 
net sediment deposition for each cell from the flow parameters. 

[0101] Step 507 involves recording the net sediment deposition for each cell. The 
record of net sediment deposition is updated within the gridding method used for the 
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simulation, typically using the Hydro-Grids gridding method described in co-pending 
U.S. Provisional Patent Application No. 60/584617. Since the horizontal lines in 
Hydro-Grids represent time lines, the cell heights are selected automatically to 
efficiently represent both thick and thin layers. Gridding methods other than Hydro- 
5 Grids may be utilized. However, these methods commonly use cells of fixed vertical 
dimension which inadequately sample very thin features and over-sample large 
homogeneous features, consuming computer memory and consuming manipulation of 
large data structures, which is time intensive. 

[0102] Step 508. Steps 505-507 are repeated until a pre-determined stopping 
10 criterion is reached. The stopping criterion is typically either the completion of a 
specified number of time steps or a specified duration of simulated time. 
Alternatively, the stopping criterion could be that the simulated deposit exceeds a 
specified height. Possible pre-determined stopping criteria include but are not limited 
to total time, thickness of the sediments, gradient of the flow, accommodation 
15 volume, and any combination thereof. 

[0103] When the inlet property values in step 504 are not known with certainty, 
an iterative process can be used where steps 504-508 are repeated, adjusting the inlet 
properties with each repetition until at least one additional stratigraphic surface 
identified in both a three-dimensional seismic data volume and in the simulated 
20 deposit are substantially similar. This additional stratigraphic surface can be 
identified by the techniques described in step 501. The similarity between the 
stratigraphic surfaces may typically pertain to their similar height above the bottom 
surface and similar lateral extent. 

[0104J An additional application of simulated deposits is to provide additional 
25 constraints and input for conventional geologic modeling methods. Such 
conventional methods include geocellular methods using conventional or multi-point 
geostatistics, object based methods, where the distribution of object shapes and sizes 
is extracted from the simulation, and spectral component methods where the high 
frequency information is derived from the simulation. The application of deposits 
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simulated according to the inventive method to the calibration of other geologic 
modeling techniques is also within the scope of this invention. 

[01051 An additional application of simulated deposits is in identifying the 
physical processes associated with creation of a particular feature of interest observed 
5 in a natural deposit. This application involves identifying a similar feature in the 
simulated deposit, and observing the flow properties associated with formation of the 
similar feature in the simulated deposit. 

EXAMPLE 

[01061 A hypothetical example of an embodiment of the inventive method will 
10 now be discussed. Figure 6 shows a three-dimensional initial bottom surface 61 
synthesized for this example, as specified in step 501 of Figure 5. Typically, the 
initial bottom surface 61 would be an interpreted stratigraphic surface obtained from 
seismic data. The initial bottom surface could also be obtained from well log data, 
core data, and/or outcrop studies. 

15 [01071 After the initial bottom surface is determined, the inlet location is 
determined (step 502 of Figure 5), a grid is created (step 503), the inlet properties are 
estimated (step 504), and the process of repeating steps 505-507 begins. Figure 7 is a 
map view of the simulation in progress, specifically showing the magnitude of the 
flow velocity (grayscale). Lighter shades indicates faster flow velocities and darker 

20 shades indicate slower flow velocities. The inlet 73, as chosen in step 502, is at the 
top center of the figure. The grid 71, as chosen in step 502, are shown superimposed 
on the image. The grid divides the simulation region into cells so that the model 
equations can be solved numerically. 

[01081 During the simulation, the net sedimentation in each cell is recorded at 
25 each time step, as specified in step 507 of figure 5. Figure 8 shows the cross-section 
of the sedimentation record once the simulation completes. Iso-time lines 81 are 
shown in black, and the median grain-size is shown in grayscale 83. Optionally, the 
simulation could be repeated after altering the inlet flow properties and/or initial 
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bottom surface until the final simulated deposit best resembles an actual deposit of 



interest. 
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Table 1 



Parameter 



A 
B 



Cd 

C 

Qo 
Cto 



Dfloc 



D 



D m 



£m 
F r 



Gi 



ho 

TO 

R 

u, 



Numeric 
Reference 
In Figures 



Definition Of Parameter 



8 



Constant, typically 1.3xl0' 7 



Empirical correction factor, typically 30. 
Half-width of the inlet. 



Bottom drag coefficient for flow, typically 0.004. 
Total volume fraction of the sediments in the sediment-laden 
flow layer. This is the sum of Ci over all grain-size bins. 



Volume fraction of the sediments of the ith grain-size bin in 

the sediment-laden flow layer. 

Volume fraction of the sediments of the ith grain-size bin in 
the sediment-laden flow layer at the inlet. 
Total volume fraction of the sediments in the sediment-laden 
flow layer at the inlet. This is the sum of Qo over all grain- 
size bins i. 



The grain-size below which the effect of flocculation 

becomes significant. 

Grain-size 



Effective grain-diameter of the ith grain-size bin 
Geometric mean grain-size suspended in the flow 



grain- 

Grain-diameter at the 10th percentile. 
Grain-diameter at the 50th percentile. 
Effective grain-size. 
Maximum value of the dimensionless erosi on rate 
Erosion rate for the ith grain-size bin 
Froude number 



27 



-— ■ — 2 

Gravitational constant, 9.8 m/s 
Volume fraction of the sediments of the ith grain-size bin in 
the surface layer of the deposit. Sum of G\ over all i is 1. 
Height of the sediment-laden flow layer. 



26 



Height of the sediment-laden flow layer at the inlet, constant 
across the width of the in let. 
Turbulent kinetic energy 

Coefficient that relates the depth averaged sediment 
concentration to the concentration at the bottom 
Ratio of the density difference between the sediment and 
water to the density of water 
Richardson number 
Particle Reynolds number. 

x- velocity component of the sediment-laden flow layer. 
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Pw 



y-velocity component of the sediment-laden flow layer. 
Velocity of the sediment-laden flow layer at the inlet. The 
velocity is assumed to be constant laterally across the width 
of the inlet and v ertically thr ough the height of the flow. 
Settling velocity 

Settling velocity o f grains of diameter Dj in still water 
Shear velocity 
Hindering function 
Effective gra in-siz e in phi scale 
Scale Factor 

Another turbulent model para meter. 
Stratification coefficient 
Turbulent model parameter. 



A constant assumed typically to be 2. 
The rate of water detrainment per unit area. 
The rate of water entrainment per unit area. 
Dimensionless height above bed. 



Model paramet er used to adjust the degree of stratification 
Bed elevation 



Von Karman constant 
A coefficient related to erosion function defined in equation 

27; 

A function related to hindered settling defined in equation 
211. 

Mean grain-size measured in phi scale. 



Kinematic viscosity of water. May be assumed to be 10" 
m 2 /s. 



Turbulent diffusion coefficient. 
Vertical profile for s ediment concentration. 
Water density. 



Sediment density. 



Standard deviation of -sizes measured in phi scale. 
Divergence of the concentration due to turbulent diffusion 
for sediments of grain-size bin i. 



Divergence of water mass as part of the turbidity flow due to 
turbulent diffusion. 

Divergence of turbulent energy due to the turbulent 
diffusion. 



Effective shear stress in x direction caused by turbulent 
diffusion. 

Effective shear stress in y direction caused by turbulent 
diffusion. 

Depositional porosity, typically 0.5 
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CLAIMS 

We claim: 

5 1 . A method for simulating formation of a sedimentary deposit comprising: 

(a) determining properties of a flow which formed a sedimentary deposit 
by solving a two-dimensional time-dependent map view system of equations for a 
sediment laden water flow, where the two-dimensional time-dependent map view 
system of equations relates flow properties including at least flow momentum, flow 

10 height, suspended sediment concentration, and entrainment of overlying clear water; 

(b) calculating net sediment deposition for at least one location using the 
flow properties from step (a); 

(c) recording the net sediment deposition from step (b) as a function of 
time for the at least one location. 

15 2. The method of claim 1 wherein the two-dimensional time-dependent map 
view system of equations is derived by assuming that vertical profiles of sediment 
concentration have similar shape at all locations and by further assuming that the 
vertical profiles of flow velocity have similar shape at all locations. 

3. The method of claim 2 wherein the vertical profiles are assumed to have a 
20 center of mass less than half the flow height. 

4. The method of claim 1 wherein the net sediment deposition is calculated as a 
difference between a deposition rate in still water and a flow velocity-dependent 
erosion rate. 

5. The method of claim 4 wherein the deposition rate in still water is determined 
25 from a still water particle settling velocity and a total sediment concentration. 

6. The method of claim 1 wherein the recording of the net sediment deposition 
from step (c) is made relative to an initial bottom topography, the initial bottom 
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topography is an interpreted surface determined from data selected from one of 
seismic data, well log data, core data, and any combination thereof. 

7. The method of claim 1 wherein the suspended sediment concentration and net 
sediment deposition is modeled for at least two ranges of grain-diameter. 

5 8. The method of claim 1 wherein the net sediment deposition is calculated using 
information from a simulation of creation and transport of turbulent kinetic energy. 

9. The method of claim 1 wherein at least one flow property is determined using 
information from a simulation of turbulent diffusion. 

1 0. The method of claim 1 wherein recording the net sediment deposition includes 
10 recording properties selected from net deposition rate for each grain-size bin, 

deposited sediment volume for each grain-size bin, grain-size distribution, bedding 
surfaces, sediment type, and any combination thereof. 

11. A method for simulating the formation of sedimentary deposits comprising: 
(a) determining an initial bottom topography; 

15 (b) determining the location of an inlet supplying a depositing flow above 

the initial bottom topography; 

(c) generating a grid dividing at least a portion of a sedimentary basin into 

cells; 

(d) estimating flow properties at the inlet; 

20 (e) determining properties of a flow forming a sedimentary deposit in at 

least one grid cell by solving a two-dimensional time-dependent map view system of 
equations for a current time step wherein the two-dimensional time-dependent map 
view system of equations relates the flow properties including at least the flow 
momentum, flow height, suspended sediment concentration, and rate of entrainment 

25 of overlying clear water, 
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(f) calculating net sediment, deposition within each cell for the current 
time step using the flow properties from step (e); 

(g) recording the net sediment deposition within each cell for the current 
time step; and 

5 (h) repeating steps (e)-(g) until a pre-determined stopping criterion is 

reached. 

12. The method of claim 11 wherein the inlet flow properties in step (d) are 
determined by an iterative process where the inlet flow properties are adjusted and 
steps (e) - (h) are repeated until the net sediment deposition is substantially similar to 

10 a measured deposit. 

13. The method of claim 12 further comprising iteratively adjusting the bottom 
topography at step (d). 

14. The method of claim 12 wherein the similarity between the simulated and 
measured deposit is determined by: 

15 identifying at least one surface in the measured deposit above the initial 

bottom topography; 

identifying at least one surface in the simulated deposit corresponding to the at 
least one surface in the measured deposit; and 

comparing the measured and simulated deposits using a criteria selected from 
20 one of shape, size, separation between corresponding surfaces of the measured and 
simulated deposits, and any combination thereof. 

15. The method of claim 14 wherein surfaces of the measured deposit are 
determined from data selected from one of well logs, cores, outcrop studies, seismic 
data, and any combination thereof. 
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16. The method of claim 1 1 further comprising using the net sediment deposition 
to generate geostatistical parameters of rock properties and further comprising using 
the geostatistical parameters in geostatistical geologic modeling techniques. 

17. The method of claim 11 wherein the net sediment deposition is used to 
5 determine shapes and sizes of geologic objects for use in object-based geologic 

■ 

modeling. 

18. The method of claim 11 wherein the net sediment deposition is used to 
interpret processes that formed a geologic feature of interest, the interpretation of the 
geologic feature of interest comprising: 

10 identifying a feature in the simulated deposit which is similar to the geologic 

feature of interest; 

observing the flow properties associated with the formation of the similar 
feature in the simulated deposit. 

19. The method of claim 1 1 wherein an HLL solver is used to calculate the flow 
15 properties. 

20. The method of claim 1 1 wherein the flow properties further comprise an 
erosion rate, wherein the erosion rate is applied in each of the cells and is determined 
by a method selected from one of selection by the user prior to the simulation, 
automatic selection based on the local flow properties in the cell, and any combination 

20 thereof. 

21. The method of claim 11 wherein the two-dimensional time-dependent map 
view system of equations is derived by assuming that vertical profiles of sediment 
concentration have a similar shape at all locations and by further assuming that the 
vertical profiles of flow velocity have a similar shape at all locations. 

25 22. The method of claim 1 1 wherein each of the vertical profiles is assumed to 
have a center of mass less than half the flow height. 
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23. The method of claim 11 wherein the net sediment deposition is calculated 
from a difference between a deposition rate in still water and a flow velocity- 
dependent erosion rate. 

24. The method of claim 23 wherein the deposition rate in still water is 
5 determined from a still water particle settling velocity and a total sediment 

concentration. 

25. The method of claim 11 wherein the initial bottom topography is an 
interpreted surface determined from using one of seismic data, well log data, core 
data, and a combination thereof. 

10 26. The method of claim 1 1 wherein the suspended sediment concentration and 
net sediment deposition is modeled for at least two ranges of grain-diameter. 

27. The method of claim 1 1 wherein the net sediment deposition is calculated 
using information from simulating creation and transport of turbulent kinetic energy. 

28. The method of claim 1 1 wherein at least one of the flow properties is 
1 5 determined using information from simulation of turbulent diffusion. 

29. The method of claim 1 1 wherein the pre-determined stopping criterion is 
selected from one of total time, deposit thickness, gradient of flow, accommodation 
volume, and any combination thereof. 

30. The method of claim 1 1 comprising recording properties of the net sediment 
20 deposition, wherein the recorded properties of the net sediment deposition are selected 

from one of a net deposition rate in each grain-size bin, deposited sediment volume in 
each grain-size bin, grain-size distribution, bedding surfaces, sediment type, and any 
combination thereof. 

31. The method of claim 11 wherein the flow properties comprise a rate of 
25 detrainment that characterizes a reduction of the flow height due to setting of the 

sediment from the depositing flow. 
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401 



Solving a 2-D time-dependent map view system 
of equations for a deep-water turbidity current; 




Calculating the net sediment deposition at each 
map view location using the flow properties 

from the map view; 



403 



Recording the time-variability of the net 

deposition. 
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Determining the initial bottom topography; 



Determining the location of the inlet supplying 

the depositing flow; 




Generating a grid dividing the region of 
simulation into cells; 




Estimating the flow properties at the inlet; 




Updating at least the flow momentum, flow 
height, suspended sediment concentration, and 
clear water entrainment rate in each cell; 




Calculating the net sediment deposition within 
each cell at the current time step; 




Recording the net sediment deposition within 
each cell at the current time step; 



508 



Repeating steps 505-508 until a stopping 

criterion is reached. 
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